home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / randist / sphere.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  3.4 KB  |  120 lines
  1. /* randist/sphere.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_rng.h>
  23. #include <gsl/gsl_randist.h>
  24.  
  25. void
  26. gsl_ran_dir_2d (const gsl_rng * r, double *x, double *y)
  27. {
  28.   /* This method avoids trig, but it does take an average of 8/pi =
  29.    * 2.55 calls to the RNG, instead of one for the direct
  30.    * trigonometric method.  */
  31.  
  32.   double u, v, s;
  33.   do
  34.     {
  35.       u = -1 + 2 * gsl_rng_uniform (r);
  36.       v = -1 + 2 * gsl_rng_uniform (r);
  37.       s = u * u + v * v;
  38.     }
  39.   while (s > 1.0 || s == 0.0);
  40.  
  41.   /* This is the Von Neumann trick. See Knuth, v2, 3rd ed, p140
  42.    * (exercise 23).  Note, no sin, cos, or sqrt !  */
  43.  
  44.   *x = (u * u - v * v) / s;
  45.   *y = 2 * u * v / s;
  46.  
  47.   /* Here is the more straightforward approach, 
  48.    *     s = sqrt (s);  *x = u / s;  *y = v / s;
  49.    * It has fewer total operations, but one of them is a sqrt */
  50. }
  51.  
  52. void
  53. gsl_ran_dir_2d_trig_method (const gsl_rng * r, double *x, double *y)
  54. {
  55.   /* This is the obvious solution... */
  56.   /* It ain't clever, but since sin/cos are often hardware accelerated,
  57.    * it can be faster -- it is on my home Pentium -- than von Neumann's
  58.    * solution, or slower -- as it is on my Sun Sparc 20 at work
  59.    */
  60.   double t = 6.2831853071795864 * gsl_rng_uniform (r);        /* 2*PI */
  61.   *x = cos (t);
  62.   *y = sin (t);
  63. }
  64.  
  65. void
  66. gsl_ran_dir_3d (const gsl_rng * r, double *x, double *y, double *z)
  67. {
  68.   double s, a;
  69.  
  70.   /* This is a variant of the algorithm for computing a random point
  71.    * on the unit sphere; the algorithm is suggested in Knuth, v2,
  72.    * 3rd ed, p136; and attributed to Robert E Knop, CACM, 13 (1970),
  73.    * 326.
  74.    */
  75.  
  76.   /* Begin with the polar method for getting x,y inside a unit circle
  77.    */
  78.   do
  79.     {
  80.       *x = -1 + 2 * gsl_rng_uniform (r);
  81.       *y = -1 + 2 * gsl_rng_uniform (r);
  82.       s = (*x) * (*x) + (*y) * (*y);
  83.     }
  84.   while (s > 1.0 || s == 0.0);
  85.  
  86.   *z = -1 + 2 * s;        /* z uniformly distributed from -1 to 1 */
  87.   a = 2 * sqrt (1 - s);        /* factor to adjust x,y so that x^2+y^2
  88.                  * is equal to 1-z^2 */
  89.   *x *= a;
  90.   *y *= a;
  91. }
  92.  
  93. void
  94. gsl_ran_dir_nd (const gsl_rng * r, size_t n, double *x)
  95. {
  96.   double d;
  97.   size_t i;
  98.   /* See Knuth, v2, 3rd ed, p135-136.  The method is attributed to
  99.    * G. W. Brown, in Modern Mathematics for the Engineer (1956).
  100.    * The idea is that gaussians G(x) have the property that
  101.    * G(x)G(y)G(z)G(...) is radially symmetric, a function only
  102.    * r = sqrt(x^2+y^2+...)
  103.    */
  104.   d = 0;
  105.   do
  106.     {
  107.       for (i = 0; i < n; ++i)
  108.     {
  109.       x[i] = gsl_ran_gaussian (r, 1.0);
  110.       d += x[i] * x[i];
  111.     }
  112.     }
  113.   while (d == 0);
  114.   d = sqrt (d);
  115.   for (i = 0; i < n; ++i)
  116.     {
  117.       x[i] /= d;
  118.     }
  119. }
  120.